
debugDict = {
    'ModSchema': ['QAM','QAM','PSK','PSK'],
    'pulseShaping': ['rrc','square','rrc','square'],
    'txAmp': [40,40,40,40],
    'rxAmp': [20,20,20,20],
    'BER': [0.495,0.59,0.45,0.5545]
}


def plotBaseBandData(data,symbol,rxSymbol=[]):
    colors = np.zeros(len(data))
    if len(rxSymbol)>0:
        figure, axis = ploot.subplots(1, 2)
        for i in range(len(rxSymbol)):
            colors[i] = rxSymbol[i]/4
        axis[1].scatter(data.real,data.imag,c=colors)
        for i in range(len(symbol)):
            colors[i] = symbol[i]/4
        axis[0].scatter(data.real,data.imag,c=colors)
        ploot.show()
    else:
        for i in range(len(symbol)):
            colors[i] = symbol[i]/4
        ploot.scatter(data.real,data.imag,c=colors)
def Convolve(data1,data2):
    ans = np.zeros(len(data1),dtype=type(data1[0]))
    for i in range(0,len(data1)):
        for j in range(0,len(data2)):
            if i-j >=0:
                ans[i] = ans[i] + data1[i-j]*data2[j]
    return ans
def getSNRvsBER(snrMin,snrMax,step):
    snr = np.linspace(snrMin,snrMax,step)
    ber = np.zeros(len(snr))
    for i in range(len(snr)):
        rxAmpDB = 20
        pbRx = np.convolve(pb,ir,mode = 'same')
        pbRx = comm.awgn(pbRx,snr[i])
        pbRx = pbRx * (10**(rxAmpDB//20))
        bbRx = comm.downconvert(pbRx, int(samplingRate/symbolRate),carrierFreq,samplingRate)
        rxSymbols = comm.demodulate(bbRx,psk)
        ber[i] = comm.ber(comm.sym2bi(rxSymbols,m),comm.sym2bi(symbols,m))
    fig, ax = ploot.subplots(1, 1, figsize=(4, 4), layout='constrained')
    ax.set_xlabel('SNR')
    ax.set_ylabel('BER')
    ax.set_ylim(0,1)
    ax.set_xscale('log')
    ax.plot(snr,ber,'x',ls='-')
    return
def getBeamPattern(dictofDegree):
    ptrn = np.zeros((37,2))
    for i in range(len(ptrn)):
        ptrn[i][0] = -180 + i*10
        try:
            ptrn[i][1] = dictofDegree[-180 + i*10]
        except:
            ptrn[i][1] = 0
    return ptrn
def plotDiffData(rx_range_min,rx_range_max,samples,opt=True):
    rx_range = np.linspace(rx_range_min,rx_range_max,samples)
    berarr = np.zeros(len(rx_range))
    indx = 0
    for i in rx_range:
        env = pm.create_env2d(frequency=27000,tx_depth=60,rx_depth=5,depth=i,rx_range=100,bottom_absorption=0.1)
        arrivals = pm.compute_arrivals(env)
        ir = pm.arrivals_to_impulse_response(arrivals, fs=100000)
        pbRx = np.convolve(pb,abs(ir)[0:int(0.1*100000)])
        pbRx = pbRx[:len(pb)]
        #pbRx = pbRx + noise(100000,len(pb),3.6,0.1)
        #pbRx = generateUAN(pbRx,samplingRate,3.6,20)
        #pbRx = comm.awgn(pbRx,20)
        pbRx = pbRx * (10**(rxAmpDB//20))
        bbRx = comm.downconvert(pbRx, int(samplingRate/symbolRate),carrierFreq,samplingRate)
        if opt:
            rxSymbols = comm.demodulate(bbRx,getDemodConst(env))
        else:
            rxSymbols = comm.demodulate(bbRx,psk)
        print("Range: " + str(i))
        berarr[indx] = comm.ber(rxSymbols,symbols)
        print("BER: " + str(berarr[indx]))
        indx = indx + 1
        plotBaseBandData(bbRx,symbols,rxSymbol=rxSymbols)
    print("Average BER: "  + str(np.sum(berarr)/len(berarr)))
def generateUAN(signal,fs,snr,w):
    #Water noise spectrum
    xf = fftfreq(len(signal),1/fs)
    nff = np.zeros(len(xf))
    nff[1:len(xf)//2] = np.sqrt(50 + 7.5*np.sqrt(w) + 20 * np.log10(xf[1:len(xf)//2]/1000) - 40* np.log10(xf[1:len(xf)//2]/1000 + 0.4))
    nff[len(xf)//2:-1] = np.sqrt(50 + 7.5*np.sqrt(w) + 20 * np.log10(-1*xf[len(xf)//2:-1]/1000) - 40* np.log10(-1*xf[len(xf)//2:-1]/1000 + 0.4))
    nff = np.power(10,nff/20)
    noise_array = np.zeros(len(signal))
    nSignal = comm.awgn(signal,snr)
    noise_array = nSignal - signal
    noise_array_fft = fft(noise_array)
    fnoise_fft = noise_array_fft * nff
    noise = ifft(fnoise_fft).real
    return signal + noise
def noise(sampRate, lenSig, w, s): 
    f = np.linspace(1, int(sampRate/2)-1, int(lenSig))/1e3 # divide by 1000 to convert Hz to kHz
    nt = 17 - 30*np.log(f) # turbulence noise
    ns = 40 + 20*(s - 0.5) + 26*np.log(f) - 60*np.log(f+0.03) #ships noise
    nw = 50 + 7.5*(w**0.5) + 20*np.log(f) - 40*np.log(f+0.04) # wind noise
    nth = -15 + 20*np.log(f) # thermal noise
    tot_noise = nt + ns + nw + nth # total noise in db per Hz
    temp = np.power(10,tot_noise/20)
    phase = 2*np.pi*np.random.randn(lenSig)
    mag = np.sqrt(temp)
    FFT = mag * np.exp(1j*phase)
    noise = np.real(np.fft.ifft(FFT))
    return noise
def getDemodConst(env):
    arrivals = pm.compute_arrivals(env)
    ir = pm.arrivals_to_impulse_response(arrivals, fs=100000)
    ir = ir[0:int(0.1*100000)]
    dataSize = 56*8
    m = 2
    n = 1
    samplingRate = 100000
    carrierFreq = 27000
    bitRate = 80
    symbolRate = bitRate//n
    txAmpDB = 40
    rrcTaps = 6
    #data Modulation
    data =comm.random_data(dataSize, 2)
    symbols = comm.bi2sym(data,m)
    symbols[10:14] = [0,1,0,1]
    symbols[14:18] = [0,1,0,1]
    psk = comm.psk(m)
    bb1 = comm.modulate(symbols,psk)
    pb =  comm.upconvert(bb1, int(samplingRate/symbolRate),carrierFreq,samplingRate)
    #Power Amplification
    pb = pb * (10**(txAmpDB//20))
    
    
    rxAmpDB = 0
    #Received Signal
    pbRx = np.convolve(pb,abs(ir)[0:int(0.1*100000)])
    pbRx = pbRx[:len(pb)]
    #pbRx = pbRx + noise(100000,len(pb),3.6,0.1)
    
    pbRx = pbRx * (10**(rxAmpDB//20))

    
    bbRx = comm.downconvert(pbRx, int(samplingRate/symbolRate),carrierFreq,samplingRate)
    ans = np.zeros(2,dtype='complex128')
    ans[0] = (bbRx[10] + bbRx[12] + bbRx[14] + bbRx[16])/4
    ans[1] = (bbRx[11] + bbRx[13] + bbRx[15] + bbRx[17])/4
    print("_______________________________DEMOD CONST__________________________________________________________________")
    plotBaseBandData(bbRx[10:18],[0,1,0,1,0,1,0,1])
    plotBaseBandData(bbRx,symbols)
    print("____________________________________________________________________________________________________________")
    return ans
def demodulate(signal,constellation):
    symbls = np.zeros(len(signal))
    for i in range(len(signal)):
        min_indx = -1
        min_val = 1000
        for j in range(len(constellation)):
            tmp = signal[i] - constellation[j]
            tmp = tmp.real ** 2+ tmp.imag**2
            if min_val > tmp:
                min_indx = j
                min_val = tmp
        symbls[i] = min_indx
    return symbls


import arlpy.uwapm as pm
import arlpy.comms as comm
import arlpy.plot as plt
import arlpy.signal as signal
import numpy as np
from scipy.fft import fft, fftfreq,ifft
import matplotlib.pyplot as ploot
import math
import matplotlib.axes as axes


pm.models()
rx_range = 1000
tx_depth = 40
rx_depth = 10
print(np.rad2deg(np.arctan((tx_depth-rx_depth)/rx_range)))
env = pm.create_env2d(frequency=27000,tx_depth=60,rx_depth=45,depth=1000,rx_range=500,tx_directionality=getBeamPattern({0:20,10:20,20:20,-20:20,-10:20,30:20,-30:20}),bottom_absorption=0.1)
arrivals = pm.compute_arrivals(env)
rays = pm.compute_eigenrays(env)
pm.plot_rays(rays, env=env, width=900)
ir = pm.arrivals_to_impulse_response(arrivals, fs=100000)
pm.print_env(env)
plt.plot(abs(ir)[0:int(0.1*100000)],fs = 100000)

1.7183580016554572

                name : arlpy
   bottom_absorption : 0.1
      bottom_density : 1600
    bottom_roughness : 0
   bottom_soundspeed : 1600
               depth : 1000
        depth_interp : linear
           frequency : 27000
           max_angle : 80
           min_angle : -80
              nbeams : 0
            rx_depth : 45
            rx_range : 500
          soundspeed : 1500
   soundspeed_interp : spline
             surface : None
      surface_interp : linear
            tx_depth : 60
   tx_directionality : [[-180.    0.]
                        [-170.    0.]
                        [-160.    0.]
                        [-150.    0.]
                        [-140.    0.]
                        [-130.    0.]
                        [-120.    0.]
                        [-110.    0.]
                        [-100.    0.]
                        [ -90.    0.]
                        [ -80.    0.]
                        [ -70.    0.]
                        [ -60.    0.]
                        [ -50.    0.]
                        [ -40.    0.]
                        [ -30.   20.]
                        [ -20.   20.]
                        [ -10.   20.]
                        [   0.   20.]
                        [  10.   20.]
                        [  20.   20.]
                        [  30.   20.]
                        [  40.    0.]
                        [  50.    0.]
                        [  60.    0.]
                        [  70.    0.]
                        [  80.    0.]
                        [  90.    0.]
                        [ 100.    0.]
                        [ 110.    0.]
                        [ 120.    0.]
                        [ 130.    0.]
                        [ 140.    0.]
                        [ 150.    0.]
                        [ 160.    0.]
                        [ 170.    0.]
                        [ 180.    0.]]
                type : 2D


print(len(ir))

# Number of samples in normalized_tone
N = len(abs(ir)[0:int(0.1*100000)])

yf = fft(ir[0:int(0.1*100000)])
xf = fftfreq(N, 1 / 100000)

plt.plot(xf, np.abs(yf))

110891


#data Generation and Config

dataSize = 2000
m = 2
n = 1
samplingRate = 100000
carrierFreq = 27000
bitRate = 80
symbolRate = bitRate//n
txAmpDB = 40
rrcTaps = 6
#data Modulation

symbols =comm.random_data(dataSize, 2)
#symbols = comm.bi2sym(data,m)
psk = comm.psk(m)
bb1 = comm.modulate(symbols,psk)
pb =  comm.upconvert(bb1, int(samplingRate/symbolRate),carrierFreq,samplingRate)

#Power Amplification

pb = pb * (10**(txAmpDB//20))

#Plotting

time = signal.time(pb,100000)
plt.plot(time,pb)


rxAmpDB = 0
#Received Signal
pbRx = np.convolve(pb,abs(ir)[0:int(0.1*100000)])
pbRx = pbRx[:len(pb)]
pbRx = pbRx + noise(100000,len(pb),3.6,0.1)
#pbRx = generateUAN(pbRx,samplingRate,3.6,20)
#pbRx = comm.awgn(pbRx,20)
#Power Amplification

pbRx = pbRx * (10**(rxAmpDB//20))

#plotting
time = signal.time(pbRx,100000)
plt.plot(time,noise(100000,len(pb),3.6,0.1))


#Demodulating Signal
bbRx = comm.downconvert(pbRx, int(samplingRate/symbolRate),carrierFreq,samplingRate)
rxSymbols = demodulate(bbRx,getDemodConst(env))

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________


bbRx = comm.downconvert(pbRx, int(samplingRate/symbolRate),carrierFreq,samplingRate)
rxSymbols = demodulate(bbRx,getDemodConst(env))

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________


comm.ber(comm.sym2bi(rxSymbols,m),comm.sym2bi(symbols,m))

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[11], line 1
----> 1 comm.ber(comm.sym2bi(rxSymbols,m),comm.sym2bi(symbols,m))

File ~\miniconda3\lib\site-packages\arlpy\comms.py:117, in sym2bi(x, m)
    115 x = _np.asarray(x, dtype=_np.int64)
    116 if _np.any(x < 0) or _np.any(x >= m):
--> 117     raise ValueError('Invalid data for specified m')
    118 y = _np.zeros((len(x), n), dtype=_np.int64)
    119 for i in range(n):

ValueError: Invalid data for specified m


plotBaseBandData(bb1,symbols)


plotBaseBandData(bbRx,symbols,rxSymbols)


plotDiffData(100,1500,50)

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 100.0
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 128.57142857142858
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 157.14285714285714
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 185.71428571428572
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 214.28571428571428
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 242.85714285714286
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 271.42857142857144
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 300.0
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 328.57142857142856
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 357.14285714285717
BER: 0.0


plotDiffData(50,150,20,opt=False)

Range: 50.0
BER: 0.0

Range: 55.26315789473684
BER: 0.0

Range: 60.526315789473685
BER: 0.875

Range: 65.78947368421052
BER: 0.0

Range: 71.05263157894737
BER: 0.133

Range: 76.31578947368422
BER: 1.0

Range: 81.57894736842105
BER: 0.0

Range: 86.84210526315789
BER: 0.0

Range: 92.10526315789474
BER: 0.0

Range: 97.36842105263159
BER: 0.0


plotDiffData(50,150,100)

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 50.0
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 51.01010101010101
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 52.02020202020202
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 53.03030303030303
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 54.04040404040404
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 55.05050505050505
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 56.06060606060606
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 57.07070707070707
BER: 0.5815

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 58.08080808080808
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 59.09090909090909
BER: 0.0


58.16326530612245
59.183673469387756


env = pm.create_env2d(frequency=27000,tx_depth=50,rx_depth=1,depth=100,rx_range=58.16326530612245,bottom_absorption=0.1)
arrivals = pm.compute_arrivals(env)
ir = pm.arrivals_to_impulse_response(arrivals, fs=100000)
plt.plot(abs(ir),fs=samplingRate)
rays = pm.compute_eigenrays(env)
pm.plot_rays(rays, env=env, width=900)


env = pm.create_env2d(frequency=27000,tx_depth=50,rx_depth=1,depth=100,rx_range=59.183673469387756,bottom_absorption=0.1)
arrivals = pm.compute_arrivals(env)
ir = pm.arrivals_to_impulse_response(arrivals, fs=100000)
plt.plot(abs(ir),fs=samplingRate)
rays = pm.compute_eigenrays(env)
pm.plot_rays(rays, env=env, width=900)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[85], line 1
----> 1 comm.ber(comm.sym2bi(rxSymbols,m),comm.sym2bi(symbols,m))

File ~\miniconda3\lib\site-packages\arlpy\comms.py:417, in ber(x, y, m)
    415     raise ValueError('Invalid data for specified m')
    416 if m == 2:
--> 417     return ser(x, y)
    418 if m > _MAX_M:
    419     raise ValueError('m > %d not supported' % (_MAX_M))

File ~\miniconda3\lib\site-packages\arlpy\comms.py:397, in ser(x, y)
    395 y = _np.asarray(y, dtype=_np.int64)
    396 n = _np.product(_np.shape(x))
--> 397 e = _np.count_nonzero(x^y)
    398 return float(e)/n

ValueError: operands could not be broadcast together with shapes (5000000,) (2000,)

0.025


#Getting Equalizer coefficients
trainSymbols = [0,1,2,3,2,1,0,2,3,1,2,0,3,1,3,0,0,1,2,3,2,1,0,2,3,1,2,0,3,1,3,0]
psk = comm.psk(4)
bb1 = comm.modulate(trainSymbols,psk)
#Upsampling
bb1 =  comm.upconvert(bb1, int(samplingRate/symbolRate),0,samplingRate)

#Power Amplification
bb1 = bb1 * (10**(txAmpDB//20))

#Channel Modelling
bbChn = Convolve(bb1,ir)
#modulating
pbChn =  comm.upconvert(bbChn,1,carrierFreq,samplingRate)
pbChn = comm.awgn(pbChn,20)
pb = comm.upconvert(bb1,1,carrierFreq,samplingRate)
#Plotting
time = signal.time(pbChn,samplingRate)
plt.plot(time,pbChn)

time = signal.time(pb,samplingRate)
plt.plot(time,pb)

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[18], line 12
      9 bb1 = bb1 * (10**(txAmpDB//20))
     11 #Channel Modelling
---> 12 bbChn = Convolve(bb1,ir)
     13 #modulating
     14 pbChn =  comm.upconvert(bbChn,1,carrierFreq,samplingRate)

Cell In[2], line 10, in Convolve(data1, data2)
      8 for i in range(0,len(data1)):
      9     for j in range(0,len(data2)):
---> 10         if i-j >=0:
     11             ans[i] = ans[i] + data1[i-j]*data2[j]
     12 return ans

KeyboardInterrupt:


#Demodulating Signal
bbRx = comm.downconvert(pbChn, int(samplingRate/symbolRate),carrierFreq,samplingRate)
bbTrueRx = comm.downconvert(pb, int(samplingRate/symbolRate),carrierFreq,samplingRate)
rxSymbols = comm.demodulate(bbTrueRx,psk)


def getEqCoeff(bbhat,bb,ntap,stepfactor=0.01,convergance=0.001):
    bbRxReal = bbhat.real
    bbRxImag = bbhat.imag
    bbRxTrueReal = bb.real
    bbRxTrueImag = bb.imag
    eqCoeffReal = np.zeros(len(bbRxReal))
    eqCoeffImag = np.zeros(len(bbRxImag))
    
    while()


ir1 = np.zeros(100000,dtype=type(ir[0]))


ir1[0] =complex(1,1)
ir1[50000] = complex(0.01,0.05)
ir1[10000] = complex(0.2,0.05)
ir1[15000] = complex(0.01,0.05)
ir1[7000] = complex(0.3,0.5)
ir1[60000] = complex(0.01,0.05)


ir1[40000] = complex(1,1)


print(len(ir1))

# Number of samples in normalized_tone
N = len(ir1)

yf = fft(ir1)
xf = fftfreq(N, 1 / 100000)

plt.plot(xf, np.abs(yf))

100000


plt.plot(abs(ir1),fs = 100000)


x = complex(1,1)

x

(1+1j)


type(x)

complex


1/40

0.025


len(beampattern[0])

2


x = {0:20

20


plt.plot(abs(ir),fs = 100000)


plt.plot(abs(ir)[0:10000],fs = 100000)


abs(ir)

[4.92515005e-02 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
 0.00000000e+00 2.02418723e-05]


k = abs(ir)[0:5000]


for i in range(5000):
    if k[i] != abs(ir)[i]:
        print(i)


plt.plot(k,fs = 100000)


!jupyter nbconvert --to html equalizer.ipynb

[NbConvertApp] Converting notebook equalizer.ipynb to html
[NbConvertApp] Writing 12068339 bytes to equalizer.html